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^ ABSTRACT 

3 ' 

We derive a fast way for measuring primordial non-Gaussianity in a nearly full-sky 
map of the cosmic microwave background. We find a cubic combination of sky maps 
combining bispectrum configurations to capture a quadratic term in primordial fluctu- 
ations. Our method takes only iV 3 / 2 operations rather than iV 5 / 2 of the bispectrum 
' analysis (1000 times faster for / = 512), retaining the same sensitivity. A key com- 

ponent is a map of underlying primordial fluctuations, which can be more sensitive to 
■ the primordial non-Gaussianity than a temperature map. We also derive a fast and 

accurate statistic for measuring non-Gaussian signals from foreground point sources. 
. The statistic is 10 6 times faster than the full bispectrum analysis, and can be used to 

^ e— | estimate contamination from the sources. Our algorithm has been successfully applied 

to the Wilkinson Microwave Anisotropy Probe sky maps by Komatsu et al. (2003). 

6 

> 

53 . 1. INTRODUCTION 

Measurement of statistical properties of the cosmic microwave background (CMB) is a direct 
test of inflation. Simple models of inflation predict Gaussian primordial fluctuations generated by 
ground-state quantum fluctuations of a scalar field (Guth Sz Pi 1982; Starobinsky 1982; Hawking 
1982; Bardeen et al. 1983; Mukhanov et al. 1992). 

Non-linearity (Salopek & Bond 1990, 1991; Gangui et al. 1994; Gupta et al. 2002), inter- 
actions of scalar fields (Allen et al. 1987; Falk et al. 1993), or deviation from the ground state 
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(Lesgourgues et al. 1997; Martin et al. 2000; Contaldi et al. 1999; Gangui et al. 2002) can gen- 
erate weak non-Gaussianity. Acquaviva et al. (2003) and Maldacena (2003) have calculated the 
2nd-order perturbations during inflation to show that simple inflation based on a slowly rolling 
scalar field cannot generate detectable non-Gaussianity with the Wilkinson Microwave Anisotropy 
Probe (WMAP) or the Planck experiments; thus, any detection of the primordial non-Gaussianity 
strongly constrains inflation models, and sheds light on physics in the early universe. For ex- 
ample, isocurvature fluctuations (Linde k Mukhanov 1997; Peebles 1997; Bucher k Zhu 1997), 
features in a scalar-field potential (Kofman et al. 1991; Wang k Kamionkowski 2000), or a "curva- 
ton" mechanism (by which late-time decay of a scalar field generates curvature perturbations from 
isocurvature fluctuations (Mollerach 1990; Lyth k Wands 2002)) can generate stronger, potentially 
detectable, non-Gaussianity. The second-order gravity also contributes to non-Gaussianity (Luo k 
Schramm 1993; Munshi et al. 1995; Pyne k Carroll 1996; Mollerach k Matarrese 1997), and there 
is a possibility that we can detect it with the Planck experiment. 

Many of the non-Gaussian models are written as a parametrized form for the curvature per- 
turbations <E>, 

= $ L (x) + f NL H(x) - (<&?»)] , (1) 

where <&l are Gaussian linear perturbations. Note that <E> = §h in Bardeen (1980). A similar 
model may also apply to isocurvature perturbations, S (Bartolo et al. 2002). This ansatz provides 
a model quantifying the amplitude of the primordial non-Gaussianity. An exact prediction of this 
model for the CMB bispectrum exists (Komatsu k Spergel 2001), while an approximate one for the 
trispectrum (Okamoto k Hu 2002). Non-linearity in inflation gives Inl ~ ©(lO" 1 ), the second- 
order gravity gives Jnl ^ C(l)> an d isocurvature fluctuations, features, or curvatons can give 
fNL 3> 1 depending on models. 

The bispectrum measured by COBE (Komatsu et al. 2002) and MAXIMA (Santos et al. 2003) 
experiments found |/jvl| ^ 10 3 (68%). Using the methods described in this paper, the WMAP 
data (Bennett et al. 2003b) have improved the constraint significantly to obtain — 58 < /nl < 134 
(95%) (Komatsu et al. 2003). While the trispectrum measured on the COBE data has shown no 
evidence for cosmological non-Gaussianity (Komatsu 2001; Kunz et al. 2001), no quantitative limit 
on fNL has been obtained. The trispectrum can be as sensitive to f^L as the bispectrum (Okamoto 
k Hu 2002); however, we need more accurate predictions using the full radiation transfer function. 
The trispectrum has not yet been measured on the WMAP data. The deficit in C 2 {6) on 9 > 60° 
(Spergel et al. 2003) might be a sign of a significant trispectrum on large angular scales. 

Measuring f^L from nearly full-sky experiments is challenging. The bispectrum analysis re- 
quires N 5 / 2 operations (iV 3//2 for computing three Vs and N for averaging over the sky) where N 
is the number of pixels (N ~ 3 x 10 6 for WMAP, 5 x 10 7 for Planck). Even though an efficient 
algorithm exists, the trispectrum still requires iV 3 (Komatsu 2001). 

Although we measure the individual triangle configurations of the bispectrum (or quadrilateral 
configurations of the trispectrum) at first, we eventually combine all of them to constrain model 
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parameters such as /nl, as the signal-to- noise per configuration is nearly zero. This may sound 
inefficient. Measuring all configurations is enormously time consuming. Is there any statistic which 
already combines all the configurations optimally, and fast to compute? Yes, and finding it is the 
main subject of this paper. A physical justification for our methodology is as follows. A model like 
equation (1) generates non-Gaussianity in real space, and central-limit theorem makes the Fourier 
modes nearly Gaussian; thus, real-space statistics should be more sensitive. On the other hand, 
real-space statistics are weighted sum of Fourier-space statistics, which are often easier to predict. 
Therefore, we need to understand the shape of Fourier-space statistics to find sensitive real-space 
statistics, and for this purpose it is useful to have a specific, physically motivated non-Gaussian 
model, compute Fourier statistics, and find optimal real-space statistics. 



2. RECONSTRUCTING PRIMORDIAL FLUCTUATIONS FROM 
TEMPERATURE ANISOTROPY 

We begin with the primordial curvature perturbations <!> (a;) and isocurvature perturbations 
S (x). If we can reconstruct these primordial fluctuations from observed CMB anisotropy, AT(ra)/T, 
then we can improve sensitivity to primordial non-Gaussianity. We find that the harmonic coeffi- 
cients of CMB anisotropy, a\ m = T _1 f d 2 nAT(n)Y^(n), are related to the primordial fluctuations 

as 

ai m = hj r 2 dr [$ im (r)a^(r) + S, ro (r)aj ao (r)] + n lm , (2) 

where &i m (r) and Si m (r) are the harmonic coefficients of the fluctuations at a given comoving 
distance, r = \x\. A beam function 6; and the harmonic coefficients of noise n\ m represent instru- 
mental effects. Since noise can be spatially inhomogeneous, the noise covariance matrix (ni m n^ m ,) 
can be non-diagonal; however, we approximate it with ~ (TQ6u'6 mm '. We thus assume the "mildly 
inhomogeneous" noise for which this approximation holds. The function ai(r) is defined by 

a i(r) = IJ VdkgnWjiikr), (3) 

where gTi(k) is the radiation transfer function of either adiabatic (adi) or isocurvature (iso) per- 
turbations. 

Next, assuming that 3> (x) dominates, we try to reconstruct $ (x) from the observed AT(n). A 
linear filter, C/(r), which reconstructs the underlying field, can be obtained by minimizing variance 
of difference between the filtered field Oi{r)ai m and the underlying field $/ m (r). By evaluating 

\O l {r)a lm -$> lm {r)\ 2 ) = U, (4) 



dOi{r) 

one obtains a solution for the filter as 



Odr) = ^, (5) 
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where the function (3i(r) is given by 

A(r) = ^ y tfdkPWffnWMkr), (6) 

and -P(fc) is the power spectrum of <£. Of course, one can replace $ with S when S dominates. We 
use a calligraphic letter for a quantity that includes effects of bi and noise such that Ci = Cft? + cTq, 
where Q is the theoretical power spectrum that uses the same cosmological model as gTi(k). 

Finally, we transform the filtered field C;(r)a; m back to pixel space to obtain an Wiener- 
filtered, reconstructed map of 3>(r, n) or S(r,n). We have assumed that there is no correlation 
between $ and S. We will return to study the case of non-zero correlation later (§ 3.2). 

Figure 1 shows Oi{r) as a function of I and r for (a) an adiabatic SCDM (fi m = 1), (b) an 
adiabatic ACDM (Q m = 0.3), (c) an isocurvature SCDM, and (d) an isocurvature ACDM. Note 
that we plot 0;(r) = (3i(r)/Ci where C\ does not include beam smearing or noise. While we have 
used P(k) oc k~ 3 for both adiabatic and isocurvature modes, specific choice of P(k) does not 
affect Oi very much as P(k) in (3[ in the numerator approximately cancels out P (k) in C\ in the 
denominator. On large angular scales (smaller V) the Sachs-Wolfe (SW) effect makes Oi —3 for 
adiabatic modes and —5/2 for isocurvature modes of the SCDM models (Sachs & Wolfe 1967). 
For the ACDM models the late-time decay of gravitational potential makes this limit different. 
Adiabatic and isocurvature modes are out of phases in /. 

The figure shows that 0\ changes the sign of the fluctuations as a function of scales. This 
indicates that acoustic physics at the last scattering surface modulates fluctuations so that hot 
spots in the primordial fluctuations can be cold spots in CMB for example. Therefore, the shape of 
0\ "deconvolves" the sign change, recovering the phases of fluctuations. This is an intuitive reason 
why our cubic statistic derived below [Eq. (9)] works, and it proves more advantageous to measure 
primordial non-Gaussianity on a filtered map than on a temperature map. 

This property should be compared to that of real-space statistics measured on a temperature 
map. We have shown in Komatsu & Spergel (2001) that the skewness of a temperature map is 
much less sensitive to the primordial non-Gaussianity than the bispectrum, exactly because of the 
cancellation effect from the acoustic oscillations. The skewness of a filtered map, on the other hand, 
has a larger signal-to-noise ratio, and a more optimal statistic like our cubic statistic derived below 
(§ 3) can be constructed. Other real-space statistics such as Minkowski functionals or peak-peak 
correlations may also be more sensitive to the primordial non-Gaussianity, when measured on the 
filtered maps; we are investigating these possibilities. 

Unfortunately, as gTi oscillates, our reconstruction of $ or S from a temperature map alone 
is not perfect. While 0\ reconstructs the primordial fluctuations very well on large scales via the 
Sachs- Wolfe effect, Oi ~ on intermediate scales (I ~ 50 for adiabatic and I ~ 100 for isocurvature), 
indicating loss of information on the phases of the underlying fluctuations. Then, toward smaller 
scales, we recover information, lose information, and so on. Exact scales at which 0\ ~ depend 
on r and cosmology. A good news is that a high signal-to-noise map of the CMB polarization 
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Fig. 1. — Wiener niters for the primordial fluctuations applied to a CMB sky map, 0;(r) = A(r)/C; 
[Eq. (5)]. We plot (a) O x for an adiabatic SCDM (Q m = 1, Q A = 0, fi b = 0.05, h = 0.5), (b) for 
an adiabatic ACDM (fi m = 0.3, = 0.7, b = 0.04, h = 0.7), (c) for an isocurvature SCDM, and 
(d) for an isocurvature ACDM. The filters are plotted at five conformal distances r = c(tq — r) as 
explained in the bottom-right panel. Here r is the conformal time (to at the present). The SCDM 
models have ctq = 11.84 Gpc and CTd ec = 0.235 Gpc, while the ACDM models ctq = 13.89 Gpc 
and CTdec = 0.277 Gpc, where r^ ec is the photon decoupling epoch. 
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anisotropy will enable us to overcome the loss of information, as the polarization transfer function 
is out of phases in I compared to the temperature transfer function, filling up information at which 
0\ ~ 0. In the other words, the polarization anisotropy has finite information about the phases of 
the primordial perturbations, when the temperature anisotropy has zero information. 



3. FAST CUBIC STATISTICS 

3.1. Primordial Non-Gaussianity 

Using two functions introduced in the previous section, we construct a cubic statistic optimal 
for the primordial non-Gaussianity. We apply filters to a/ m , and then transform the filtered cty m 's 
to obtain two maps, A and B, given by 

A(r,n) = ^2^i aim Y lm (n), (7) 

lm 

B( r ,n) = ^2^P^a lm Y lm (n). (8) 

Ira Ll 

The latter map, B(r,n), is exactly the Oj-filtered map, an Wiener-filtered map of the underlying 
primordial fluctuations. We then form a cubic statistic given by 

S W im = 4* J r 2 dr j ^A(r,n)B 2 (r,n), (9) 

where the angular average is done on full sky regardless of sky cut. We find that S pr i m reduces 
exactly to 

gobs gprim 

prim ~ L> c h c h c h ' 1 j 

h<h<h 1 2 3 

where 

&hhh = B hhh b h b h b h, (11) 

and B°bf 2l3 is the observed bispectrum with the effect of 6/ corrected while the theoretical 

one for /jvl = 1 (Komatsu & Spergel 2001), 



where 



= 2I hl2l3 J r 2 dr [(3 h (r)Pl 2 (r)a l3 (r) + l3 (r)P h (r)a h (r) + p l2 (r)P h (r)a h (r)] , (12) 



W = V (2 " + " ( \ +1)(2h + 1) (o o o)- <13> 

It is straightforward to derive equation (10) from (9) using equation (12). 

The denominator of equation (10) is the variance of Bf^f 2l in the limit of weak non-Gaussianity 
(say \f NL \ < 10 3 ) when all Vs are different: (B 2 ihh ) = C h C l2 C h A hhh , where A hhh is 6 for 
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h = h = h, 2 for h = h ^ h etc., and 1 otherwise. The bispectrum configurations are thus 
summed up nearly optimally with the approximate inverse-variance weights, provided that Aj 1 j 2 j 3 
is approximated with ~ f . The least-square fit of Bf^™ to Bf^f 2l can be performed to yield 

Sprim^fNL £ ( 14 ) 



This equation gives an estimate of Jml directly from S, 



prim,' 



The most time consuming part is the back-and-forth harmonic transformation necessary for 
pre-filtering [Eqs. (7) and (8)], taking N 3 ^ 2 operations times the number of sampling points of r, of 
order 100, for evaluating the integral [Eq. (9)] . This is much faster than the full bispectrum analysis 
which takes N 5 / 2 , enabling us to perform a more detailed analysis of the data in a reasonable 
amount of computational time. For example, measurements of all bispectrum configurations up 
to Imax = 512 take 8 hours to compute on 16 processors of an SGI Origin 300; thus, even only 
100 Monte Carlo simulations take 1 month to carry out. On the other hand, S pr i m takes only 30 
seconds to compute, 1000 times faster. When we measure /nl for l max = 1024, we speed up by a 
factor of 4000: 11 days for the bispectrum vs 4 minutes for S pr i m . We can do 1000 simulations for 
lmax = 1024 in 3 days. 



3.2. Mixed Fluctuations 

The 0/-filtered map, B, is an Wiener-filtered map of primordial curvature or isocurvature 
perturbations; however, this is correct only when correlations between the two components are 
negligible. On the other hand, multi-field inflation models (Langlois 1999; Gordon et al. 2001) and 
curvaton models (Lyth & Wands 2002; Moroi & Takahashi 2001) naturally predict correlations. The 
current CMB data are consistent with, but do not require, a mixture of the correlated fluctuations 
(Amendola et al. 2002; Trotta et al. 2001; Peiris et al. 2003). In this case, the Wiener filter for the 
primordial fluctuations [Eq. (5)] needs to be modified such that Oi(r) = (3i(r)bi/Ci — > f3i(r)bi/Ci, 
where 

M d \r) = \ j\ 2 dk[p^k)g a T f{k)+P c {k)g^{k)\3l{kr), 

Pt so (r) = \ j k 2 dk [P s (k)gff(k) + P c (k)g a T f(k)\ j t (kr), 

for curvature (adi) and isocurvature (iso) perturbations, respectively. Here P<j> is the primordial 
power spectrum of curvature perturbations, Ps of isocurvature perturbations, and Pq of cross 
correlations. 

As for measuring non-Gaussianity from the correlated fluctuations, we use equation (1) as a 
model of 3>- and 5-field non-Gaussianity to parameterize them with fff[ and /^|£, respectively. We 
then form a cubic statistic similar to S pr i m [Eq. (9)], using A(r,n) and a new filtered map B(r,n) 
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which uses 0i(r). We have two cubic combinations: A a ^B^ di for measuring fff[ and Ai so Bf so for 
f^l, each of which comprises four terms including one Pj (or Pj), one P^, and two PjPc's (or 
PsPc's). In other words, the correlated contribution makes the total number of terms contributing 
to the non-Gaussianity four times more than the uncorrelated-fluctuation models (see (Bartolo 
et al. 2002) for more generic cases). 



3.3. Point Source Non-Gaussianity 

Next, we show that the filtering method is also useful for measuring foreground non-Gaussianity 
arising from extragalactic point sources. The residual point sources left unsubtracted in a map can 
seriously contaminate both the power spectrum and the bispectrum. We can, on the other hand, 
use multi-band observations as well as external template maps of dust, free-free, and synchrotron 
emission, to remove diffuse Galactic foreground (Bennett et al. 2003a). The radio sources with 
known positions can be safely masked. 

The filtered map for the point sources is 

D(n) = J2waimYi m (n). (15) 

This filtered map was actually used for detecting point sources in the WMAP maps (Bennett et al. 
2003a). Using D(n), the cubic statistic is derived as 



/j2 ~ o teobs t2src 

h<h<h 1 2 3 



Here, Bf^ is the point-source bispectrum for unit white-noise bispectrum, = h^i^- (b S rc = 

1 in Komatsu & Spergel (2001).) When covariance between Pf^*™ an d Bf^ is negligible as is the 
case for WMAP and Planck (Komatsu & Spergel 2001), we find 

We omit the covariance only for simplicity; including it is simple (Komatsu & Spergel 2001; Komatsu 
et al. 2002). Aga in S src measures b src much faster than the full bispectrum analysis, constraining 
effects of residual point sources on CMB sky maps. Since S src does not contain the extra integral 
over r, it is even 100 times faster to compute than S pr i m . This statistic is particularly useful because 
it is sometimes difficult to tell how much of C\ is due to point sources. Komatsu et al. (2003) have 
used S src ( 

i.e., b src ) to measure C\ due to the unsubtracted point sources. 
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3.4. Incomplete Sky Coverage 



Finally, we show how to incorporate incomplete sky coverage and pixel weights into our statis- 
tics. Suppose that we weight a sky map by M(n) to measure the harmonic coefficients, 



A full-sky ai m is related to through the coupling matrix Miy mm i = f d 2 hM{h)Y lm (h)Yii m i{h) 
by a im = Yli'm' a i'm'^Wmm' ■ m this case the observed bispectrum is biased by a factor of 
J d 2 hM 3 (n) / (4tt); thus, we need to divide S pr i m and S src by this factor. If only sky cut is consid- 
ered, then this factor is a fraction of the sky covered by observations (Komatsu et al. 2002). 

We have carried out extensive Monte Carlo simulations of non-Gaussian sky maps computed 
with equation (2). Appendix A of Komatsu et al. (2003) describes the simulations in detail. We 
find that S pr i m reproduces input Jnls accurately both on full sky and incomplete sky with modest 
Galactic cut and inhomogeneous noise expected for WMAP, i.e., the statistic is unbiased. We 
cannot however make a sky cut very large, e.g., more than 50% of the sky, as for which the covariance 
matrix of &\ x i 2 i z is no longer diagonal. The cubic statistic does not include the off-diagonal terms 
of the covariance matrix [see Eq. (10)]; however, it works fine for WMAP sky maps for which we 
can use more than 75% of the sky. Also, we have found that equation (17) correctly estimates b src 
using simulated realizations of point sources (see Appendix B of Komatsu et al. (2003)). As for 
uncertainty in the cubic statistics, Figure 2 shows that errors of /jvl from S pr i m and b src from S src 
are as small as those from the full bispectrum analysis (see descriptions in Appendix). 



Using the method described in this paper, we can measure non-Gaussian fluctuations in a 
nearly full-sky CMB map much faster than the full bispectrum analysis without loss of sensitivity 
(see Appendix). Our fast statistics allow us to carry out extensive Monte Carlo simulations charac- 
terizing effects of realistic noise properties of experiments, sky cut, foreground sources, and so on. 
A reconstructed map of the primordial fluctuations, which plays a key role in our method, poten- 
tially gives other real-space statistics more sensitivity to primordial non-Gaussianity. As we have 
shown our method can be applied to the primordial non-Gaussianity arising from inflation, gravity, 
or correlated isocurvature fluctuations, as well as the foreground non-Gaussianity from radio point 
sources, all of which can be important sources of non-Gaussian fluctuations on the forthcoming 
CMB sky maps. 




(18) 



4. 



CONCLUSION 
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A. 



PERFORMANCE OF CUBIC STATISTICS 



We have extensively tested our cubic statistics using Monte Carlo simulations of CMB sky 
with realistic noise properties and various galaxy masks. More specifically, we have used WMAP 
1-year noise properties in V band (Bennett et al. 2003b) and straight masks in Galactic coordinates 
with \b\ cut = 0, 10, 20, 30, 40, 50, 60, 70, and 80°. (I.e., -\b\ cut < b < \b\ cut has been masked.) 
Figure 2 shows uncertainty in /nl and b src obtained from 300 Gaussian simulations using the cubic 
statistics (diamonds) for different galaxy masks. The solid lines show the minimum variance which 



would be expected for the full bispectrum analysis. These lines have been computed by y F~ / f s ky, 
where Fij is the Fisher matrix (Komatsu &: Spergel 2001) and f s k y is a fraction of sky surviving 
the mask. We find that the cubic statistics perform remarkably well, giving errors as small as 
the full bispectrum analysis for all masks, for simulations of the CMB signal only and CMB with 
homogeneous noise. (The homogeneous noise has r.m.s. noise which matches the average noise in 
V band.) However, the statistics perform a bit worse in the presence of inhomogeneous noise - this 
is because we are using a simple uniform weighting (M(n) = for the masked pixels and M(n) = 1 
otherwise) rather than the optimal C~ l weighting. Since noise in some pixels are much larger than 
average, our statistics are affected by those noisier pixels. It is not straightforward to implement 
C^ 1 weighting (which is non-diagonal) in our cubic statistics, but we continue to investigate a way 
to take into account inhomogeneous noise in our method. 
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Fig. 2. — Performance of cubic statistics. The left panels and right panels show errors of /jvl 
and bgrc, respectively, which are obtained from 300 Gaussian simulations. Each point has been 
computed for a given straight sky cut with |6| CMi . From the right to left, \b\ cu t = 0, 10, 20, 30, 
40, 50, 60, 70, and 80°. The solid lines show the minimum variance which would be obtained by 
the full bispectrum analysis. The top, middle, and bottom panels show simulations of the CMB 
signal only, CMB plus homogeneous noise, and CMB plus inhomogeneous noise, respectively. Noise 
properties assume WMAP 1-year data in V band. 
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